library(tidyverse)
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.1 ──
## ✔ ggplot2 3.4.0 ✔ purrr 0.3.4
## ✔ tibble 3.1.7 ✔ dplyr 1.0.9
## ✔ tidyr 1.2.0 ✔ stringr 1.4.0
## ✔ readr 2.1.2 ✔ forcats 0.5.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
library(here)
## here() starts at /Users/jasonmoggridge/Documents/CUBE_UOttawa
library(ggiraph)
library(patchwork)
source(here('R/report_plots.R'))
##
## Attaching package: 'lubridate'
## The following objects are masked from 'package:base':
##
## date, intersect, setdiff, union
swabs <- read_rds(here('data/cube.rds')) |>
filter(str_detect(site, '^UO_')) |>
glimpse()
## Rows: 721
## Columns: 39
## $ date_swab <date> 2021-09-21, 2021-09-21, 2021-09-21, 2021-09-21…
## $ city <fct> Ottawa, Ottawa, Ottawa, Ottawa, Ottawa, Ottawa,…
## $ site <fct> UO_90U, UO_90U, UO_90U, UO_DMS, UO_DMS, UO_DMS,…
## $ floor <chr> "Level 1", "Level 1", "Level 1", "Level 1", "Le…
## $ location <chr> "Site 1: Entrance Foyer, Just Before The Securi…
## $ replicate <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ lot <chr> "XG024", "XG024", "XG024", "XG024", "XG024", "X…
## $ barcode <chr> "10349", "2921", "10310", "14573", "11841", "10…
## $ pcr_positive <fct> Negative, Negative, Negative, Negative, Negativ…
## $ worker <fct> Aaron, Alex, Rees, Partha, Naila, Keaton, Rache…
## $ lab <fct> Ottawa, Ottawa, Ottawa, Ottawa, Ottawa, Ottawa,…
## $ type <fct> University, University, University, University,…
## $ outbreak <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ susp_outbreak <fct> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ outbreak_flag <fct> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ comments <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ room_demographics <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ cases <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ case_dates <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ case_comment <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ date_pcr <date> 2021-09-24, 2021-09-24, 2021-09-24, 2021-09-24…
## $ pcr_result <chr> "not detected", "not detected", "not detected",…
## $ negative_control <lgl> FALSE, FALSE, TRUE, FALSE, FALSE, TRUE, FALSE, …
## $ spoiled <lgl> FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE…
## $ co2 <dbl> 538, 519, 538, 582, 528, 582, 730, 687, 730, 54…
## $ date_shipped <date> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
## $ pcr_result2 <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ pcr_result3 <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ location_within_room <fct> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ comment <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ case_comments <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ redo <lgl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ swab_type <chr> "flock", "flock", "flock", "flock", "flock", "f…
## $ case_role <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ pcr_ct <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ floor_location <chr> "Level 1: Site 1: Entrance Foyer, Just Before T…
## $ sample_location <fct> "Uo_90u: Level 1: Site 1: Entrance Foyer, Just …
## $ replicate_barcode <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ replicate_swab_1_dates <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
swabs |> count(site, date_swab)
## # A tibble: 316 × 3
## site date_swab n
## <fct> <date> <int>
## 1 UO_90U 2021-09-21 3
## 2 UO_90U 2021-09-23 3
## 3 UO_90U 2021-09-28 3
## 4 UO_90U 2021-09-30 3
## 5 UO_90U 2021-10-05 2
## 6 UO_90U 2021-10-07 2
## 7 UO_90U 2021-10-12 2
## 8 UO_90U 2021-10-14 2
## 9 UO_90U 2021-10-19 2
## 10 UO_90U 2021-10-21 3
## # … with 306 more rows
# waste water data
raw_ww <-
readxl::read_excel(here::here('data/wastewater.xlsx')) |>
janitor::clean_names()
uottawa_ww <-
raw_ww |>
filter(!is.na(sample_date)) |>
mutate(
source =
source_name |>
str_remove('^Data - uOttawa - ') |>
str_remove('.xlsx$'),
site = case_when(
source == 'uo_aa' ~ 'AA',
source == 'uo_fa' ~ 'FA',
source == 'uo_ft' ~ 'FT',
source == 'uo_na' ~ 'NA',
source == 'uo_rgn' ~ 'RGN',
source == 'uo_st' ~ 'ST',
source == 'uo_tbt' ~ 'TBT',
TRUE ~ 'not matched'
)
)
# ottawa wastewater data: summary
summarise_wastewater <- function(){
read_rds(here('data/ww_ottawa.rds')) |>
select(sample_date, starts_with('cov_')) |>
pivot_longer(contains('cov_')) |>
mutate(target = str_extract(name, 'cov_n.'),
stat = str_extract(name, 'mean|sd'),
) |>
select(-name) |>
pivot_wider(names_from = stat, values_from = value) |>
mutate(.lower = mean - sd, .upper = mean + sd)
}
wifi <-
read_rds(here('data/uo_wifi.rds')) |>
filter(date >= min(swabs$date_swab), date <= max(swabs$date_swab))
This plot shows the counts of positive (red) and negative (yellow) samples collected at each facility over time (x-axis). Samples that could not be tested are shown in navy. Only flocked swabs are shown. (Other sponge swabs were collected on 2022-04-28 were for comparing flocked and sponge swabs.)
plot_timeseries_summary(
swabs |>
filter(
swab_type %in% 'flock',
!negative_control,
!is.na(negative_control)
),
height_svg = 3.3,
width_svg = 8
)
This plot shows the counts of positive and negative samples collected at each facility over time.
plot_timeseries_locations(
swabs |>
filter(swab_type == 'flock'),
height_svg = 6,
width_svg = 8
)
uottawa_ww |>
ggplot(aes(sample_date, viral_copies_per_copies_pep_avg, color = site)) +
geom_line(alpha = 0.5) +
geom_point(alpha = 0.5) +
labs(x = 'Date', y = 'SARS-CoV-2 / PPMoV')
## Warning: Removed 16 rows containing missing values (`geom_line()`).
## Warning: Removed 80 rows containing missing values (`geom_point()`).
rm(raw_ww, uottawa_ww)
This section contains results from modeling covid cases at UO using swab results as a predictor.
Derek and Jason used mixed effects logistic regression to model the occurrence of cases (y/n) based on swab results for the previous week (the proportion of positive swabs). The model has a random intercept for sites.
I switched this model to a Bayesian mixed-effects model. In the mixed-effects generalized linear model, MRT/90U are problem-sites where the random intercept parameter causes the model fit to be overspecified (model matrix becomes singular).
Our model formula is
cases ~ positives[lag 1 week] + (1|site).
The model is fit as follows:
swab_model <-
blme::bglmer(
cases_binary ~ detection_lag_1week + (1 | site),
data = uo_sites,
family = binomial
)
swab_model <- read_rds(here('model/uo_swab_model.rds'))
These plots show the swab results, cases, and predictions by the
current model.
# uo cases data - long table
cases_long <-
read_rds(here('./data/latest_uo_cases.rds')) |>
filter(!is.na(positive_result),
positive_result > lubridate::as_date('2021-09-14')) |>
select(case, role, confirmed_by_test, positive_result, c(`90U`:TBT)) |>
mutate(
`Campus Total` = TRUE,
# 3 days transmissibility
transmission_start = positive_result - days(3)
) |>
pivot_longer(`90U`:`Campus Total`, names_to = 'site') |>
filter(value) |>
select(-value) |>
mutate(site = fct_relevel(site, "Campus Total", after = 0L))
uo_pred <- read_rds(here('model/uo_pred.rds'))
# create a 3-panel fig for each of 6 buildings.
map(
# site names for filtering
.x = swabs$site |> unique() |> str_remove('UO_'),
.f = ~plot_uo_pred_swab_case(
swabs = swabs,
cases = cases_long,
preds = uo_pred,
site_select = .x
)
) |> wrap_plots(guides = 'collect', ncol = 2)
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
These tables show the model coefficients and statistics.
broom.mixed::glance(swab_model) |>
mutate_if(is.numeric, round, 2) |>
kableExtra::kable(caption = 'Model statistics',
format = "html") |>
kableExtra::kable_styling(bootstrap_options = "striped",
full_width = F, position = "left")
## Loading required package: blme
## Loading required package: lme4
## Loading required package: Matrix
##
## Attaching package: 'Matrix'
## The following objects are masked from 'package:tidyr':
##
## expand, pack, unpack
| nobs | sigma | logLik | AIC | BIC | deviance | df.residual |
|---|---|---|---|---|---|---|
| 84 | 1 | -24.52 | 55.04 | 62.33 | 42.83 | 81 |
broom.mixed::tidy(swab_model) |>
mutate(p.value = as.character(p.value)) |>
mutate_if(is.numeric, round, 3) |>
mutate(p.value = as.numeric(p.value)) |>
kableExtra::kable(caption = 'Model parameters',
format = "html") |>
kableExtra::kable_styling(bootstrap_options = "striped",
full_width = F, position = "left")
| effect | group | term | estimate | std.error | statistic | p.value |
|---|---|---|---|---|---|---|
| fixed | NA | (Intercept) | -3.301 | 0.804 | -4.107 | 0.0000401 |
| fixed | NA | detection_lag1w | 7.114 | 2.112 | 3.369 | 0.0007549 |
| ran_pars | site | sd__(Intercept) | 1.067 | NA | NA | NA |
# tidy(swab_model, effects = 'ran_vals')
sjPlot::tab_model(
swab_model,
show.re.var = TRUE,
pred.labels = c("(Intercept)", "Swabs @ t-1week"),
dv.labels = "Relation between UO cases (y/n) and proportion\n of positive swabs the previous week"
)
| Relation between UO cases (y/n) and proportion of positive swabs the previous week | |||
|---|---|---|---|
| Predictors | Odds Ratios | CI | p |
| (Intercept) | 0.04 | 0.01 – 0.18 | <0.001 |
| Swabs @ t-1week | 1229.31 | 19.59 – 77128.28 | 0.001 |
| Random Effects | |||
| σ2 | 3.29 | ||
| τ00 site | 1.14 | ||
| ICC | 0.26 | ||
| N site | 6 | ||
| Observations | 84 | ||
| Marginal R2 / Conditional R2 | 0.378 / 0.538 | ||
This plot shows the modelling data as points (detection lagged 1 week and cases at a given site- y/n), as well as how the probability of future cases varies by the previous weeks detection level and site, according to our model (curves).
read_rds(here('fig/plt_model_resp_curves.rds')) +
labs(subtitle = 'Points show case/swab data and curves show model probability values')
## Warning: Removed 42 rows containing missing values (`geom_point()`).
This plot shows the odds ratios for the intercepts of the model (an intercept for each site).
# random effect plot: site intercepts; axes are flipped
sjPlot::plot_model(swab_model, type = 're') +
theme_light() +
labs(title = "Random effects", y = 'Odds Ratio')
This panel shows the cases counts at UO over the course of our sampling period. The case data shown represents the days on which a positive test was reported (black rug lines) and the presumed start of transmissiblity for each case (red lines).
# generated from file R/plot_UO_case_load.R
read_rds(here('fig/plt_cases_swabs.rds'))
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## Warning: Removed 1 rows containing non-finite values (`stat_bin()`).
## Warning: Removed 4 rows containing missing values (`geom_bar()`).
This panel shows linked data from swab results, CO2
readings, and wifi traffic (number of users @ peak, daily) during our
study period. Unfortunately, we do not have wifi data for 90U.
source(here('R/uott_presentation.R'))
uo_multivariate & theme(text = element_text(size = 13))
## Warning: Removed 1 rows containing missing values (`geom_point()`).
rm(uo_multivariate)
plot_UO_wifi_ts(wifi, swabs)
This panel shows a time-series of the daily peak number of wifi users
at UO facilities. Sampling days are highlighted in blue.